TUTORIAL 6 : Exercise 1

Author:Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)

Reproducing U(r) with EE & WL/RE methods

In the limit of infinite dilution, i.e. when the solute concentration (or number density) is negligible, the potential of mean force approaches the underlying effective interaction for a pair of solute particles.

This implies that if only two particles are present in the simulation box,

\[{\rm PMF}=-\ln g(r)\ \to\ \beta\sum_k U_k(r)\]

In this exercise we will see that it is actually the case.

Exercise 1.1: Using the EE method

First, we will consider a pair of LJ particles in vacuum.

Change to directory tutorial6-1/fed_2ljp_ee/, where you will find the input files for this exercise.

Check the cutoff, energy units and particle diameter in FIELD file:

Lennard-Jones - 2 particles
CUTOFF 12.0
UNITS K
NCONFIGS 1
atoms 1
Lj core  1.0  0.00
MOLTYPES 1
test
MAXATOM  2
vdw 1
Lj core  Lj core  lj 1.000000 3.000
FINISH
CLOSE

Check the FED input in CONTROL file:

# FED block is started by 'use fed <flavor> <stride>' (must be in 'use' section on top)
use fed generic 1

# FED <method> <damping> <scaling> <update> [<smooth_itr> <smooth_1st> <smooth_last> <weight>]
fed method EE   0.5       1.125     500000    #3  10  250  #0.333

# FED order [param] <type> <nbins> <first> <last> <power>
fed order parameter com2    250     2.5     27.5   1
com1 mol 1 atom 1   # use COM1 for a group of molecules or atoms within molecules
com2 mol 1 atom 2   # use COM2 for a group of molecules or atoms within molecules
com  sampling correction 1 # entropy on a spherical surface (on the fly)
#com  sampling correction 2 # entropy in a spherical bin volume (a posteriori)
fed order parameter done

# close FED block started by 'use fed' directive
fed done

NOTE: Since at every MC step one of the FED-involved particles is moved, the stride for FED attempts can only be ‘1’ (it will default to that otherwise)!

We start an EE iteration with 500,000 MC steps per EE update cycle initially. After two first cycles the update stride will be doubled every time the FED bias is updated (and stored in FEDDAT.<node>_<cycle> file), i.e. we are doing a sequence of 2 x 500,000; then 1,000,000; then 2,000,000; and so on…

The EE method requires a damping factor, \(\eta\), for the energy feedback term, \(\Delta U_{bias}(R)\) (see above), which in this case is \(\eta_0=0.5\). The damping factor is scaled up after every EE iteration by the scaling coefficient as \(\eta_{k+1}=1.125\ \eta_k\) which is done with the same frequency as defined by the main update stride.

The subsection defining the order parameter is started by directive ‘fed order …’ and closed by ‘fed order [param] done’. The first (opening) record in this subsection specifies: the parameter grid number, the allowed range of variation (working window), and the power of incremental bin size (if greater than 1 in absolute value).

The specs for the two COM groups are self-evident: each contains a single atom (indeces 1 and 2) in the same ‘molecule 1’ (under name ‘test’ in FIELD and CONFIG). The third directive in the ‘fed order’ subsection specifies how the bias unfolding correction is to be accounted for. In this case we opted for the ‘on the fly’ correction that accounts for the ideal entropy on a sphere of the actual (current) radius, \(\sim 2\ln R\).

Further in CONTROL file the number of MC steps and strides for collecting stats should be commensurate with those for FED:

temperature          1.0         # Reduced temperature
nocoul all                       # No electrostatics
#distewald                       # Distributed Ewald sums in parallel runs

equilibration        0           # Equilibration period (not for FED calcs)
steps                4000000     # Total number of MC steps 2^(<iters>-1)*<update>
print                1000000     # Frequency of printing to OUTPUT
stats                 100000     # frequency of printing stats to PFILE
check                1000000     # Check rolling energies vs recalculated from scratch
sample rdfs     300 30.0 100     # Sampling grid, bin size & frequency for RDF calculation
sample coords          10000     # Frequency of storing configuration in trajectory

acceptatmmoveupdate  1000000000  # DO NOT update max atom displacement
maxatmdist           10.000      # Max atom displacement
move atoms         1 100         # Move atoms
Lj core                          # Name & type of atom to move

start simulation

Run the simulation (serial version) and check the output:

[tutorial6-1/fed_2ljp_ee]$ ls FED*
FEDDAT.000_001  FEDDAT.000_002  FEDDAT.000_003  FEDDAT.000_004

Launch gnuplot and plot the FED data, which is the second column in FEDDAT files (to quit gnuplot type ‘q’):

[tutorial6-1/fed_2ljp_ee]$ gnuplot
gnuplot> plot 'FEDDAT.000_001' u 1:2 w l t "EE iter 1", 'FEDDAT.000_002' u 1:2 t "EE iter 2" w l, \
'FEDDAT.000_003' u 1:2 w l t "EE iter 1", 'FEDDAT.000_004' u 1:2 w l t "EE iter 4"

Now let’s add the known target (LJ potential) and zoom in by providing X and Y ranges:

gnuplot> lj(arg) = 4*( (3.0/arg)**12 - (3.0/arg)**6 )
gnuplot> plot [x=2:10] [y=-1.5:1.5] 'FEDDAT.000_001' u 1:2 w l t "EE iter 1", \
'FEDDAT.000_002' u 1:2 t "EE iter 2" w l, 'FEDDAT.000_003' u 1:2 w l t "EE iter 3", \
'FEDDAT.000_004' u 1:2 w l t "EE iter 4", 'FEDDAT.000_001' u 1:(lj($1)) w l lc 'black' lt 0 t "Target U(r)"

Finally, also check the EE population histograms, i.e. probability distributions, which is the third column in FEDDAT files:

gnuplot> plot [] [y=0:12000] 'FEDDAT.000_001' u 1:3 w l t "EE iter 1", \
'FEDDAT.000_002' u 1:3 t "EE iter 2" w l, 'FEDDAT.000_003' u 1:3 w l t "EE iter 3", \
'FEDDAT.000_004' u 1:3 w l t "EE iter 4"

What might go wrong?

Questions to answer:

• Was the entire order parameter range covered appropriately?
• Are the histograms sufficiently flat?
• Can one improve of the stats by varying \(\eta_0\), the scaling-up coef, or the bias update stride?
• Would smoothing upong updating the bias help?

Refinement production run is necessary!

Let’s store away the current results:

[tutorial6-1]/fed_2ljp_ee$ mkdir equil-fed1/
[tutorial6-1]/fed_2ljp_ee$ cp CON* FIELD equil-fed1/
[tutorial6-1]/fed_2ljp_ee$ mv *.0??* equil-fed1/

Now amend CONTROL file:

# FED <method> <damping> <scaling> <update> [<smooth_itr> <smooth_1st> <smooth_last> <weight>]
fed method EE   0.7       1.125    2000000    #3  10  250  #0.333

# FED order [param] <type> <nbins> <first> <last> <power>
fed order parameter com2    250     2.5     27.5   -1
...
steps                8000000     # Total number of MC steps 2^(<iters>-1)*<update>

and re-run.

Ooops…

Due to ‘-1’ at the end of ‘fed order’ directive, we are asking for a seed file FEDDAT.000 which we forgot to provide! Of course, we meant to use the FED bias data from the latest iteration:

[tutorial6-1]/fed_2ljp_ee$ cp equil-fed1/FEDDAT.000_004 ./FEDDAT.000

restart it once again and redo the analysis of the production run.

Exercise 1.2: Using the WL method with replica-exchange in T

Next, we will apply the Wang-Landau method to the same system of a pair of LJ particles.

Change to directory tutorial6-1/fed_2ljp_wl/, where the input files are found.

We will combine the WL FED calculation with the replica-exchange (RE) parallel tempering technique so that FED(R) will be evaluated for four temperatures at once!

To this end, the following two lines in CONTROL file are added (in the ‘use’ section below the title):

use ortho                  # use orthorhombic PBC for speed-up (if the cell is cubic)
use repexch 4  -0.2  1000  # Replica-Exchange with 4 replicas, delta T = -0.2, sampling every 1000 MC steps

The FED input to use the WL method is as follows:

# FED block is started by 'use fed <flavor> <stride>' directive
use fed generic 1

# FED <method> <droplet> <scaling> <update> [<smooth_itr> <smooth_1st> <smooth_last> <weight>]
fed method WL   0.002     0.5       500000    #3  10  250  #0.333

# FED order [param] <type> <nbins> <first> <last> <power>
fed order parameter com2    250     2.5     27.5   1

Make sure that reading the FED bias from input is turned off, i.e. ‘<power>’ is set ot 1 in ‘fed order …’ directive. Also, make sure that the rest of the input in CONTROL corresponds to our initial run comprising 4,000,000 MC steps in total.

Similar to EE, the WL iteration initially starts with 500,000 MC steps per WL update cycle. The update stride is doubled in the same manner as for EE: 2 x 500,000; then 1,000,000; then 2,000,000; and so on…

The WL method requires an initial droplet value, \(\delta\) (recall the WL illustration), which in this case is \(\delta_0=0.002\). The scaling coefficient acts as follows: \(\delta_{k+1}=0.5\ \delta_k\) with the same frequency as defined by the main update stride.

NOTE: Even though there are no charges (no electrostatic interactions) in this case, one has to invoke the ‘distewald’ directive for any parallel run!

temperature          1.0         # Reduced temperature
nocoul all                       # No electrostatics
distewald                        # Distributed Ewald sums in parallel runs

equilibration        0           # Equilibration period (not for FED calcs)
steps                8000000     # Total number of MC steps 2^(<iters>-1)*<update>
print                1000000     # Frequency of printing to OUTPUT
stats                 100000     # frequency of printing stats to PFILE
check                1000000     # Check rolling energies vs recalculated from scratch
sample rdfs     300 30.0 100     # Sampling grid, bin size & frequency for RDF calculation
sample coords          10000     # Frequency of storing configuration in trajectory

acceptatmmoveupdate  1000000000  # DO NOT update max atom displacement
maxatmdist           10.000      # Max atom displacement
move atoms         1 100         # Move atoms
Lj core                          # Name & type of atom to move

start simulation

To run the replica exchange parallel job, you have to ensure that you are not running an interactive session, but instead submit the job to a parallel queue:

[tutorial6-1/fed_2ljp_wl]$ sbatch parallel.sub

Run the simulation and check the output (order can be different):

[tutorial6-1/fed_2ljp_wl]$ ls FED*
FEDDAT.000_001  FEDDAT.000_002  FEDDAT.000_003  FEDDAT.000_004  FEDDAT.000_005
FEDDAT.001_001  FEDDAT.001_002  FEDDAT.001_003  FEDDAT.001_004  FEDDAT.001_005
FEDDAT.002_001  FEDDAT.002_002  FEDDAT.002_003  FEDDAT.002_004  FEDDAT.002_005
FEDDAT.003_001  FEDDAT.003_002  FEDDAT.003_003  FEDDAT.003_004  FEDDAT.003_005

Launch gnuplot and plot the FED data:

[tutorial6-1/fed_2ljp_wl]$ gnuplot
gnuplot> plot [x=2:10] [y=-3:2] 'FEDDAT.000_005' u 1:2 w l t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:2 t "Replica 2, WL iter 5" w l, 'FEDDAT.002_005' u 1:2 w l t "Replica 3, WL iter 5", \
'FEDDAT.003_005' u 1:2 w l t "Replica 4, WL iter 5"

Check the histograms too:

[tutorial6-1/fed_2ljp_wl]$ gnuplot
gnuplot> plot [] [y=14000:18000] 'FEDDAT.000_005' u 1:3 w l t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:3 t "Replica 2, WL iter 5" w l, 'FEDDAT.002_005' u 1:3 w l t "Replica 3, WL iter 5", \
'FEDDAT.003_005' u 1:3 w l t "Replica 4, WL iter 5"

• Store away the data, amend the input and run the refinement (production) simulation on your own!

Extra exercises

• Try to run EE/RE combination

• Try to add charges (+1,-1) on the LJ particles and redo the FED calculations - what do you expect to obtain?

• Study a more physically meaningful case:

  TUTORIAL 6 : Exercise 2 - Calculating PMF for two charged colloids in ionic cloud

  TUTORIAL 6 : Exercise 3 - Umbrella sampling in subranges (windows) with the use of WHAM method